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Abstract 

In this paper wc present the parton level Monte Carlo program TeVJet, a direct im- 
plementation of the dipole subtraction method for calculating jet cross sections in NLO 
QCD. It has been written so as to allow the inclusion of new processes in as straight- 
forward a way as possible. The user must provide the usual ingredients for an NLO 
calculation and from these the process-independent parts required to make the phase 
space integrals finite in 4 dimensions are automatically generated. These integrals are 
then performed using Monte Carlo techniques. We present the results for a few example 
processes. 
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1 Introduction 



Since the equations of motion derived from the QCD Lagrangian are non-linear coupled equa- 
tions, they cannot be solved exactly. At large momentum transfers we can exploit the prop- 
erty of asymptotic freedom and use perturbation theory to expand about the partonic states 
of quarks and gluons and calculate partonic cross sections. In the case of hadron colliders 
one must then fold these with non-perturbative but universal parton distribution functions. 
One of the most powerful tests of QCD and our understanding of the strong interaction has 
been the comparison between experimental data and such perturbative calculations. 

Leading order (LO) calculations can be useful in interpreting the experimental data, and 
there exist several programs that have automated the calculation of observables to this order 
(some examples are given in Refs. [1-3]). Higher order terms in perturbation theory have a 
richer structure, indeed it is not until one calculates the 1-loop correction that we encounter 
ultraviolet divergences and need to renormalize the field theory. Leading order calculations 
can only, for some scattering processes, describe the shapes of distributions well but do not 
say much about the overall normalization. This manifests itself as a strong dependence on 
the unphysical renormalization and factorization scales. 

Thus the calculation of observables to next to leading order (NLO) accuracy is more 
involved, which, at a practical level, is due to the presence of singularities at intermediate 
steps in the calculation. This includes both the ultraviolet divergences already mentioned, 
which require the machinery of renormalization to tame; and infrared divergences, which 
are a consequence of a field theory with massless quanta. These infrared singularities pose 
a significant practical challenge in the calculation of exclusive observables, i.e. when the 
available phase space is restricted. The problem is that the cancellation occurs between two 
quantities that result from integration over phase space volumes of different dimensionality, 
which makes the numerical integration over these phase space volumes non-trivial — each of 
these integrals is divergent (although their sum is finite). 

There are basically two methods for overcoming these problems and performing a (semi-) 
numerical calculation of these phase space integrals: the phase space slicing method and 
the subtraction method. Both of these exploit the universality of the singular behaviour of 
tree level QCD amplitudes in the soft and coUinear limits. One of the advantages of such 
numerical calculations is that several different observables (jet functions) can be calculated 
in parallel. 

The earliest NLO calculation of an event shape distribution [4] introduced and used the 
subtraction method. In Ref. [5], this was developed to produce the first general-purpose 
program, designed to calculate the NLO corrections for many observables simultaneously, 
for a single process, e^e~ annihilation. Other early calculations used the phase space slicing 
method, with varying degrees of success [6]. Refs. [7,8] defined the first general method for 
calculating arbitrary infrared-safe observables in arbitrary processes, using the phase space 
slicing method. Refs. [9,10] did the same thing for the subtraction method, defining the dipole 
subtraction method, which has been the basis for the majority of numerical NLO calculations 
performed since then. Ref. [10] also succeeded in solving the analogous, but more general, 
problem with an arbitrary number of identified particles in the final state. This was also 
done for the phase-space slicing method in Ref. [11]. More recently, those seeking to define 
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a subtraction algorithm for NNLO calculations (for examples, see Refs. [12-14]) have been 
lead to define new variations on the dipole subtraction theme. The generalization to massive 
partons was done for the phase-space slicing method in Ref. [11] and the dipole subtraction 
method in Ref. [15]. 

In this publication we describe the parton level Monte Carlo program TeVJet which is 
a direct implementation of the dipole subtraction method, at present in its massless version 
without identified particles in the final state [10]. In section [5] we give a brief overview of 
the dipole subtraction method; in section [3] we describe in detail the structure of the code; 
in section [3] we describe how to use the program with a simple example main program and 
in section [5] we present some results for the test processes e^e~ — > 3 jets and ep — > 1 jet. 



2 The Dipole Subtraction Method 
2.1 Overview of the method 

In this section we introduce the basic idea of the dipole subtraction method [10] for the 
simple case of processes with no hadrons in the initial state (i.e. e^e~ annihilation processes). 
In section 12.21 we shall generalize this to include processes with hadron beams. Finally in 
sections 12. 3112. 51 we shall discuss some of the details that are useful for an understanding of 
TeVJet. 

At next-to-leading order (NLO) in perturbative QCD some observable, a, is given by 

where a^^ is given by the integration of the Born level cross section over the available phase 
space, i.e. 



da"". (2) 
The subscript signifies that this is an integral over m-body phase space so explicitly: 



^Mnj;-i- /#(™)|A^^|2Ff\ (3) 

{rn\ ' 



where !F is the flux factor. A/in represents all of the QCD independent factors, X]{m} denotes 
the sum over all m-parton configurations, Sj^}. is the Bose symmetry factor, dcj)'^'^'^ is an 
element of m-body phase space, is the matrix element in the Born approximation, 

and Fj^^ is the jet function, which defines the observable to be calculated. The jet function 
warrants a little explanation - it is the value of the observable as a function of the momenta of 
the m final state particles. If one wishes to calculate a total cross section this could be simply 
given by Fj = 1, or for a total cross section with some cuts it could be some combination of 
theta functions such that it is given by Fj = 1 for all kinematic configurations satisfying the 
cut, and -Fj = for all other configurations. For a differential cross section it will contain 
a delta function. Alternatively it could simply contain an analytical function of the final 
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state momenta. In general it will contain a combination of these possibilities but there are 
some constraints on it that we shall come back to later. Note that by definition the leading 
order contribution, CT^O, is finite and the phase space integration can be carried out in four 
dimensions. This implies that either refers to an intrinsically infrared safe process, or 
the jet function vanishes {Fj"^^ — > 0) in all soft and collinear limits {pi-pj 0). 

The next-to- leading order contribution, cr'^^*-', is given by the sum of two parts - the real 
part, which is the contribution from the real emission Feynman diagrams, and the virtual 
part, which is the interference between the one-loop and Born level amplitudes: 

J m+l J in 

^+l)|_^i^|2^j"^+l)^ j d4>^^)\M^\''Ff'\ (5) 

where is the real matrix element, and by slight abuse of notation we write the interference 
term as These are separately divergent in four dimensions although their sum is finite. 

In fact the sum of these two terms is only finite for infrared safe observables, defined as an 
observable whose value is independent of the number of soft and collinear hadrons (partons) 
in the final state. This means that we require the following property of the jet function: 

pirn+l) ^ pirn) ^ (6) 

as one approaches an m + 1-parton configuration that is kinematically degenerate with an 
m-parton configuration. 

One way to proceed is to regularize these integrals, e.g. by dimensional regularization, 
which would lead to single (1/e) and double (1/e^) poles, and then calculate the integrals 
analytically. Once the two are summed the e poles will cancel and the limit e ^ can 
be taken. However we would like to calculate the integrals numerically using Monte Carlo 
techniques, so we need another strategy. 

The behaviour of tree level QCD amplitudes in the soft and collinear limits is well known. 
The function |A4^p is a complicated function that diverges in several regions of phase space, 
whenever a gluon becomes arbitrarily soft, or when any pair of massless partons becomes 
collinear. The basic idea of the dipole subtraction algorithm is to construct a local counter 
term for the real contribution: 

NLO _ [ u^R j^^i , / , / 



a''^^= / / / da^ (7) 

J m+l J m+l Jm 

such that the first term on the RHS of Eq. ([7]) is finite. The counter term, da^, is constructed 
as the sum of several functions each of which diverges in fewer places than da^. In fact it is 
constructed, using an improved factorization, from the Born level cross section: 

da^= da^ dVdipoie, (8) 

dipoles 

where by da^ we mean an appropriate colour + spin projection of the Born level cross section; 
and the dipole term, dVdipoie, represents a 2 — > 3 splitting process. 
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One should think of these terms as being produced in a two step process: first an m- 
parton configuration with Born-Uke kinematics is produced; then one of these particles (the 
emitter) decays into two particles in the presence of another (the spectator). The dipole 
terms are process independent and universal and only depend on the quantum numbers of 
the emitter and spectator. As we shall see in the practical implementation, the construction 
of the counter term is the reverse of this process - first an m + 1-parton configuration is 
produced and then this is mapped onto an m-parton configuration. 

Of course the remaining two integrals on the RHS of Eq. ([Tj) must be combined for the 
method to work, and we shall sketch out how this works now. In Eq. ([H]) the subtraction 
cross section, da^, is evaluated at the m + 1-parton phase space point {pi, . . . ,pm+i)- It 
is constructed by choosing three of the m + 1 partons, i,j,k (summing over all possible 
such choices), and projecting their momenta, Pi,Pj,Pk onto new momenta, which we call 
Pij,Pk, and all other momenta {pi} onto {pi}- That is, formally, the tilde momenta, {p}, are 
functions of the momenta {p} and of which partons are involved in the splitting process, i.e. 
z, j and k. 

We describe ij and k as the emitter and spectator respectively. Their momenta are chosen 
to be linear combinations of Pi,Pj,Pk such that: they are on mass-shell {pfj = p| = 0); 
momentum is conserved {p^ + Pj + p^ = pfj +Pk)] and the phase space of the partons i,j, k 
can be factorized into the dipole phase space times a single parton phase space: 

d<j)^^\pi,Pj,Pk;Q) = d(l)^^\pij,pk;Q)[dpi{pij,pk)] 

= Jd^^^\pij,pk;Q)dn^^-^^dzdy. (9) 

Here Q is the total momentum of the 3-parton system, ^ is a Jacobian factor for the mapping, 
and z and y are the degrees of freedom of the {d — l)-dimensional single parton 

subspace. In particular, ^}^'^~^^ represents a set of d — 3 angles that define the direction of 
parton i in the transverse subspace: in four dimensions it reduces to the azimuthal angle (p. 
All other momenta remain unchanged. The second term on the RHS of Eq. ([7]) is given by: 



/ da' 

Jm+l 



[ #(-^)(p; Q) ^ (10) 

J dip^ks'^'^ [dpi{Pij,Pk)] 



Since the three parton phase space factorizes and the dipole term, dV^ipoie/dpi, is the only 
term that depends on the single parton phase space we can writ^ 

Jm+1 diprfes-^ #(™)(p) J [dpiiPij,Pk)] 



(12) 



^Note a slight subtlety in going from Eq. (|10|l to Eq. (I12|l : Exchanging the order of summing and integrating 
has changed the meaning of the momenta; In Eq. (I10|l . each dipole defines a different set of tilde momenta 
{p} for the same momenta {p}, whereas in Eq. pip , each dipole defines a different set of momenta {p} for the 
same tilde momenta {p}. Since the momenta are dummy variables of integration and the exact phase space 
factorization has the same form for each dipole, this is a completely transparent change. Having made it, the 
order can be re-exchanged as in Eq. p2|l . 
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where the insertion operator, I{p; e), is given by: 




The insertion operator contains all of the divergent behaviour of da and can be calculated 
analytically in d dimensions. Since the dipole term, Vdipoio and the Jacobian factor, J, are 
process independent this can be calculated once and for all. 

This can now be combined with the virtual term as it is an integral over m-parton phase 
space (see section . The NLO contribution is therefore given by: 

^NLO _ ^NLO{m+l} _|_ ^NLO{m} 

= / [da^-da^]+ f [da^ + da^0l]. (14) 

Now we have two phase space integrals, which are finite in 4 dimensions and can be calculated 
by Monte Carlo integration. 



2.2 Initial state hadrons 

When there are hadrons in the initial state, the perturbatively calculable partonic cross 
section must be folded with non-perturbative, but universal, parton distribution functions. 
That is the cross section for a reaction with two hadrons in the initial state with momenta 
pi and p2 is given by: 

Cr{Pl,P2) = yZ dr]afh/a{Va,tJ'F) / dr]b fh/b{Vb, fJ-p) 
a,b 

X Wab^iVaPl,ribP2) + crait'^{riaPi,ribP2, ^f)], (15) 

where the sum is over parton flavours; the parton distribution function (pdf), fh/aiVa, fJ-'p), 
gives the probability of finding a parton type a inside the hadron with a fraction ija of the 
momentum of the hadron at some scale /x^; o"^^ and c^J"^ are the partonic cross sections, 
which we shall define below. 

The leading order contribution is given by: 

= / -^MnY.^d<P^^^M^b\\VaPuVbP2)F;r\ (16) 

where ric and account for the average over the number of initial state spin and colour states, 
is the flux factor. A/in accounts for all QCD independent factors, M.^^ amplitude, 
in the Born approximation, for partons type a and b to give an m-parton final state that 
contributes to the observable; the sum X]{m} runs over all such final states. 

The presence of hadrons (and hence partons) with well defined momenta in the initial 
state spoils the cancellation of divergences between the real and virtual contributions. The 
parts in the real contribution that do not cancel are the initial-state collinear divergences. 



6 



i.e. when a final-state parton becomes coUinear with one of the initial-state partons. The 
problem is in fact one of double counting: this is exactly the physics that is already included 
in the pdfs; and the solution is to subtract off these collinear poles from the real term, 
which introduces a scale dependence into both the pdfs and the cross section. The NLO 
contribution is therefore given by: 

c^ah^{jPi,V2,lil)= f da^,,{pi,p2)+ [ da^^ipi,p2)+ f dag{pi,p2, f^l), (17) 

where the last term is the collinear subtraction term. Although the pole must be subtracted, 
there is a free choice of how much finite physics is also subtracted and absorbed into the 
pdf. A particular choice of this sets the factorization scheme and by changing the scheme 
we simply move physics between the pdf and the perturbative calculation. The collinear 
counter-term is given by: 

d(Tab{Pi,P2,fi'l) = ~|^ r(l - e) ^ dz2da^{zipi,Z2P2) 

1 /A^..2 



5uS{l - Z2) 

+ 5ac5{l - Zi 



e 



-^(^)^^^^(^2) + i^^sm, (18) 



where P°-^{z) are the (four-dimensional) regularized Altarelli-Parisi splitting functions and 
the kernels K'^^a,{z) specify the factorization scheme (setting Kp^{z) = defines the MS 

subtraction scheme). 

In the case of initial state partons there are additional phase space constraints due to the 
fact that there are hadrons in the initial state with well defined momenta. These constraints 
are diff'crcnt depending on whether the emitter and spectator are in the final state or initial 
state. Therefore there are different dipole contributions: 

da^= J2 '^'^^ <^ (^^dipole + d^dipole)' (19) 
dipoles 

where (iVdipoic is the dipole term already discussed for emitter and spectator both in the 
final-state, and dV^^^^^^ represents the remaining three dipole terms involving initial state 
emitters and/or spectators. 

As before we consider the analytical integration of each dipole contribution over the one- 
parton subspace leading to the divergence. For the integration of the new dipoles involving 
initial state partons the analogue of the phase space factorization is a phase space convolution. 
The specific form of the convolution depends on the phase space constraints, but in the case 
of a final state emitter, ij, and an initial state spectator, a, it has the form: 

d(t)^'^\pi,Pj\Q +Pa) = / dxd(t)^'^\pij\Q + XPa)[dpi{pij\Pa,x)\ 

Jo 

= / dxd(t)^^\pij;Q + xpa)Jdn^'^-^Uz, (20) 
Jo 
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where here J is a new Jacobian factor. This dipole contribution is given by: 



da 



B 



d^M [dpi{pij;pa,x)] 



dx / + 



da 



B 



dV°-- 



[dpi{pij;Pa,x)] 



.(21) 



The last part contains all of the dependence on the single parton subspace variables 0^*^"^^ 
and z, and is process independent. Thus this integral can be carried out once and for all in d 
dimensions. In a similar way the contributions from the other two types of dipole with initial 
state emitters and/or spectators can be calculated. There is a part of each of these terms 
that is proportional to 5(1 — x) and these terms combine with the integral of the final-state 
dipole, Vdipoie, to form the insertion operator, I(p;e). This matches the singular behaviour 
of the virtual term. 

It can be shown that all of the divergences from the remaining terms cancel with those of 
the collinear counter-term in Eq. (|18p . Thus all of the 1/e poles cancel and the finite parts 
can be collected together to give a convolution term of the form: 



dx 



<^ab {x;xpa,Pb,fJ'F) + (^ab [.X;Pa, Xpb, ^Ip) 



(22) 



where 

-I 







dxa'^l''^"'\x-xpa,Pb,t^l) = Y, [ dx I \daa'b{xPa,Pb)®{'P + 'i^r'{x)]^^^, (23) 



and analogously for a^l''^^'^''^ {x;pa, xph, /x|,). Therefore after the dipole subtraction procedure 
has been carried out the NLO contribution to the jet observable may be written in the 
following form: 



NLO{m+l} 



'ah 



NLO{m} 
ab 



+ 



dx 



a 



NLO{m} 
ab 



, -2 \ , .NLO{m} / -2 N 

{x;xpa,Pb,f^F) + (^ab [x; Pa, XPb, f^p) 



(24) 



where each of the three terms is finite in four dimensions and can be calculated by Monte 
Carlo integration. 



2.3 Factorization in the soft and collinear limits 

In order to construct the dipole subtraction terms, we need to know the behaviour of tree- level 
QCD amplitudes in the divergent regions of phase space, i.e. in the soft and collinear limits. 
In this section we recall these properties of tree-level QCD amplitudes, and in section [23] we 
shall write down the dipole terms. In order to do both of these it is necessary to introduce 
the notation we shall use [10]. 
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The tree-level matrix element for a scattering process with QCD partons in the initial- 
state and m QCD partons in the final-state i^ 

-MJ^;a::.':'"''='"'---''^'---'''"'''"'---(Pl, ...,Pm;Pa,...), (25) 

where {ci, . . . , c^, Ca, • • • }, {si, . . . , Sm, Sa, ■ ■ ■} and {pi, . . . ,pm,Pa, ■ ■ ■} label respectively the 
colour, spin and momenta of the external partons. We introduce a basis {\ci, . . . , Cm, Ca, ■ ■ ■)0 
\si, . . . , Sm, Sa, ■ ■ ■)} in helicity + colour space in such a way that 



^m;a;.'.':'"''"""''''""''"''"""(Pl> ■ ■ ■ , Pm] Pa, ■ ■ ■ ) = VMa) ■ ■ ■ 
X ((ci, . . . , 

Cm, Ca, ■ ■ ■ \ \Si, . . . 1 Sm, Sa, ■ ■ ■ \ ) \ ^, ■ ■ ■ ,'m] (1^ . . . )m, (26) 

where nda) is the number of colour states of the initial state parton a {Nc for a quark/anti- 
quark, — 1 for a gluon). That is |1, . . . , m; a, . . . )m is a vector in colour + helicity space. In 
this notation the matrix element squared summed(averaged) over final(initial)-state colours 
is given by: 

|A^m;a,...P =m;a,... {I, ■ ■ ■ , m; a, . . . \l, . . . , 171] a, . . . )m;a,...- (27) 

In order to describe colour correlations it is convenient also to define a colour charge operator, 
Tj, associated with the emission of a gluon from each parton i. If the emitted gluon has colour 
index c, the colour charge operator is 

Ti = T[\c) (28) 
and its action onto the colour space is defined by 

(ci , . . . , Cj , . . . , Cm , Ca, • • • , c|Tj 1 6i , . . . , 6j , . . . , bm ,ba, ■ ■ 

= 5c^b^ . . . T^^h^ . . . Sc^brr^^caba ■■■ , (29) 

where T^" = ifabc if the emitting particle i is a gluon and T^^ = f^^ if the emitting particle 
is a quark (in the case of an emitting anti-quark T^^ = t^^ = —t^^). In this notation 
conservation of colour is 

^^T, + T„ + ...^ |l,...,m;a,...)^ = 0. (30) 

Finally, in this notation the square of the colour-correlated matrix element needed for the 
soft limit is given by 

l^mfa,...P =m;a,... {1, . . . , m; a, . . . \T i .T j\l, . . . , m; a, . . . ) m;a,..., (31) 

where the indices I, J refer either to final-state or initial-state partons. 



^In the following we shall use ... to represent final-state partons, and a, . . . to represent initial-state 
partons. 
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2.3.1 Soft limit 



Consider the behaviour of a QCD matrix element with m + 1 final-state partons when one 
final-state gluon becomes arbitrarily soft. If the momentum of this gluon is pj, the soft limit 
can be parameterized by: 

p^^=\q^ A^O, (32) 

where is an arbitrary 4-vector. In the soft limit the matrix element squared with m + 1 
partons tends towards the matrix element squared with m partons found by removing one 
gluon, multiplied by a divergent factor 

m+i;a,...(l, . . . ,m + l;a, . . . |1, . . . ,m + l;a, . . . )m+i;a,... -^4TTfi'^''as 

Xm;a,.,.(l, . . . , m + 1; a, . . . I [3^{q)]^ J/.(9)|l, • • • , m + 1; a, . . . )m;a,...- (33) 
The eikonal current is given by 

J'^('^) = ET/^, (34) 

where the sum is over all final-state and initial-state partons. Putting these together we 
re-write the soft limit factorization as 

m+i;a,...(l, . . . ,m + l;a, . . . |1, . . . ,m + l;a, . . . )m+i;a,... ^ --^Svr^^'as 

X I]/ ^ I]E-^/m;a,...(l, . . . , m + 1; a, . . . \ ^{pj^'pi^),q' \ '^^ . . . , m + 1; a, . . . )m;a,...- 

(35) 

We point out two properties of the factorization in Eq. ()35p that are relevant for the 
construction of the dipole subtraction terms. First we note that the Born level matrix 
element squared does not factorize exactly, for there are colour correlations relating to the 
fact that an arbitrarily soft gluon cannot resolve the colour charges of individual partons, 
but is emitted coherently from (at next-to-leading order) all pairs of partons. It is one of the 
roles of the spectator parton in the dipole subtraction method to account for these colour 
correlations. Next we note that this factorization is only valid in the soft limit. If one wished 
to define a counter term for the soft limit in such a way that the factorization was valid over 
the whole of phase space some care is needed. In a 1 — > 2 splitting process, if one wishes to 
ensure that all external particles are on mass-shell, momentum cannot be conserved. This 
is the other role of the spectator in the dipole method - it absorbs some momentum so 
that momentum conservation is implemented exactly (see note in section 12.4.21 about the 
initial- initial dipole 2?"*'''). 
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2.3.2 CoUinear limit 



Now consider the limit in which two final-state QCD partons with momenta pi and pj become 
colUnear, parameterized as follows: 



z 2p.n 

1.2 ij, 



2pi.pj 



J V ± l_z2p.n 



z{i-zy 

k± 0, (36) 

where p^ is a light-like vector (p^ = 0) that gives the collinear direction; kj_ is a time-like 
vector transverse to the collinear direction {k'j_ < 0, k±.p = 0); and n'^ is an additional 
light-like vector necessary to specify how the collinear limit is approached, which satisfies 
k±.n = 0. As we approach the collinear limit the matrix element squared for m -|- 1 partons 
factorizes into the matrix element squared for m partons multiplied by a divergent splitting 
kernel: 

r- (37) 

The m parton matrix element is found by replacing partons i and j with a single parton 
ij with the addition rule (from QCD vertices) that anything-l-gluon becomes anything, and 
quark-|-antiquark becomes a gluon. In the notation of Ref. [10]: 

m+l;a,... {1, ■ ■ ■ ,m + 1; a, . . . \1, . . . ,m + 1; a, . . . )m+l;a,... -^4.-Kii^^as 

Pi-Pj 

Xm;o,...(l, ... ,m + l;a,... \P(^ij)^i{z, e)|l, . . . , m -|- 1; a, . . . )m;a,... ■ (38) 

The splitting kernel P(jj) ^(z, fcj^; e) is the unregularized d-dimensional Altarelli-Parisi 
splitting function, note that it depends not only on the collinear momentum fraction, z, 
but also on the transverse momentum k±. We note that the Born level matrix clement does 
not factorize exactly, in fact there are spin correlations arising from the fact that the splitting 
kernel is a matrix in the helicity space of the parton ij acting on the helicity indices of parton 
ij in the vectors m;a,...(l) . . . , m -|- 1; a, . . . | and |1, . . . , m -|- 1; a, . . . )m;a,...- We do not list the 
splitting functions here, but refer the reader to [10]. 

Next consider the limit that a final-state parton i becomes collinear with an initial-state 
parton a, parameterized as: 

p1 = {l-^)p'a+k1_ ^' 



^Pi-Pa 



1 — X 2pa.n 



1-x' 

k± 0, (39) 
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and the splitting process a ^ ai + i involves a transition from the initial-state parton a to 
the initial-state parton ai with the associated emission of the final-state parton i. Again the 
quantum numbers of the emitter ai are assigned according to their conservation at the tree- 
level QCD vertices: if a and i are the same type, the ai is a gluon; if a is a fermion(gluon) and 
i is a gluon(fermion) then ai will be a fermion(antifermion). The collinear limit as Pi.pa — > 
reads 

,m + l;a, ...)m+i;a,... 



Pa,(ai){z,kr,^)\^, - ■ ■ ,m + l;ai, ...)m;ai,... • (40) 

The m-parton matrix element on the RHS is found from the m + 1-parton matrix element 
by removing parton i, and replacing parton a with parton ai (note that this parton carries 
momentum xpa)- 

As with the soft limit the factorization in Eqs. (j38p and ()40p is only valid in the collinear 
limit, and some care must be taken when extrapolating away from the collinear regions of 
phase space. Finally we note that the soft and collinear limits overlap: it is possible for 
a final-state gluon to become both soft and collinear with another external parton. When 
constructing a local counter term for da^ one must ensure that one avoids double-counting 
of these singular regions. 

2.4 The dipole contributions 

In this section we shall construct the local counter term, da^, for the real term. It has the 
property that as we approach the soft and collinear limits it reproduces the factorizations of 
Eqs. (f35]) and (fSHIl . i.e. as we approach a singular area of phase space 

da^ ^ da^. (41) 

Since momentum is conserved exactly, these limits are approached smoothly thus avoiding 
double counting in the regions that are both soft and collinear. As well as requiring that da^ 
matches the singular behaviour of da^ we also require that it has no additional singularities - 
a point that we shall return to. 

In fact there are several ways of extrapolating smoothly from the singular regions of phase 
space, and we shall introduce four different counter terms to match the different phase space 
constraints. In the following we shall use a pair of indices to denote an emitter, and a single 
index to denote a spectator. These will appear as a subscript (superscript) to denote a final 
(initial) state parton. Figure [T] shows the four different types of dipole contribution that act 
as counter terms for the real matrix element squared, \Mm+i;a,...\^ in different regions of 
phase space. 



^+l;a,...(l, . . . ,m + l;a, . . . |1, 



1 1 



XPi-Pa 

Xm;ai,...(l, . . . , m + 1; ai , , 
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Figure 1: Illustration of the four dipole types we consider in this paper: final-state emit- 
ter with final-state spectator iVij^k), final-state emitter with initial-state spectator (^^i^), 
initial-state emitter with final-state spectator (P"*) and initial-state emitter with initial-state 
spectator (P''*'*). 



2.4.1 Final-state singularities 

Using the dipole factorization formula we can write the real matrix element squared as we 
approach a final-state singularity {pi-Pj — > 0) as 

m+l;a,b 

{l,...,m + l;a,6|l, ... ,m + l;a,b)m+i;a,b 

= X] ^ihk{Pl, ■ ■ ■ ,Pm+i;Pa,Pb) 

1, ■ ■ ■ ,Pm+i;Pa,Pb) +1^ij{Pl, ■ ■ ■ ,Pm+i;Pa,Pb) \ , (42) 

where we have neglected all finite terms. The first term on the RHS is a dipole contribution 
with a final-state emitter and a final-state spectator and is given by 

• • • ,Pm+i;Pa,Pb) = -ir-^ — 
2pi.pj 

Xm;a,fe(l, . . . ,ij , ■ ■ ■ ,k, . . . ,171 + 1; o, 6| ^^^^^ Vij,fc| 1, . . . ,ij , . . . ,k, . . . ,m + l;a,b)m;a,b- (43) 



The m-parton matrix element on the RHS of Eq. (I43p is obtained from the m + 1-parton 
matrix element by replacing (a) the partons i and j with a single parton ij (the emitter) and 
(b) the parton k with the parton k (the spectator). 
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The quantum numbers (except momenta) of the partons of the m-parton matrix element 
are assigned as fohows. The quantum numbers of all partons except the emitter and the 
spectator are as in the m + 1-parton matrix element. The quantum numbers of the spectator, 
k are the same as k. The quantum numbers of the emitter, ij, are obtained according to 
their conservation in the collinear splitting process ij — > i+ j, i.e. if i and j are something + 
gluon then ij is something; if i and j are quark + antiquark then ij is a gluon. 

The momenta of the emitter and spectator are defined as follows 



Pk 



1 



1 - yij,k 

where the scalar yij^k is given by 



~u u . a 

Pij =Pi +Pj 



-Z^ 

yij,k 



yij,k 



Pi-Pj 



Pi-Pj +Pj-Pk+Pk-Pi 



(44) 



(45) 



Note that the emitter and spectator are both on mass-shell and in the splitting process 
{ij,k} — > {i,j,k} momentum conservation is implemented exactly. 

Note that the matrix element squared does not factorize exactly on the RHS of Eq. ()43p , 



there are colour and spin correlations. Tjj and T^, are the colour charges of the emitter 
and spectator respectively, and Vjj ^ is an operator in the helicity space of the emitter, 
which is related to the d-dimensional Altarelli-Parisi splitting functions. As we wish only to 
outline the general method rather than go into every detail we shall not list here the splitting 
matrices for each type of dipole contribution. We shall write out these matrices for the case 
of final-state emitter and spectator, but for the other three types we refer the reader to [10]. 
These splitting functions are given in terms of the Lorentz scalars yij k, % and Zj given 

by 



Pi-Pk 



_ Pi-Pk 

Pj-Pk+Pi-Pk Pij-Pk' 
~ ^ Pi-Pk ^ Pj^k 
^ Pj-Pk+ Pi-Pk Pij-Pk 

For a fermion -|- gluon splitting we have 

2e, ^ r 2 



{s\'yq,gj,k{^i;yij,k)\s') = S-Kfl^asCp 



1 - Zi. 



(46) 



1 - Zi{l - yij^k) 



[1 + Zi) - e(l - z, 



^qiQj ,k^ss' ) 



where s and s' are the spin indices of the fermion ij in the vector | . . . , ij , 
antiquark and gluon -|- gluon splitting we have 



^ss' 

(47) 

For quark -|- 



iiv 



Pi-Pj 



-{Zipf - Zjp'^){ZiP^ - Zjp"^) 



(48) 
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1 



+ 



1 



1 



I - Zi{l- yij^k) 'i- - - yij,k) 



Pi -Pj 



(49) 



where /i and v are the spin indices of the gluon ij in the vector | . . . , , . . . ) . 

The other two terms on the RHS of Eq. (|42p are dipole contributions with a final-state 
emitter and an initial-state spectator. This is given by 



VUpi, . . . ,pm+i;Pa,Pb) 



1 1 

I .Tq, 



Xm;a,6(l, ■ ■ ■ , ij , ■ ■ ■ , m + l] a, b\ % ° V".|1, . . . , ij, . . . , m -h 1; a, 

where the Lorentz scalar Xij^a is given by 

_ Pi-Pa +Pj -Pa -Pi-Pj 



ij,a 



{Pi+Pj)-Pa 



(50) 



(51) 



All quantum numbers except the momenta of the emitter and spectator are given as in 
the case of the dipole term Pjj^fc. These are given by 



Pa = Xij,aP'i, Pii=Pi +P j - Xij,a)P^ ■ 



(52) 



For a definition of the splitting matrices V^- we refer the reader to the original paper [10]. 
2.4.2 Initial-state singularities 

Using the dipole factorization formula we can write the real matrix element squared as we 
approach an initial-state singularity {pa-Pi — > 0) as 



m+l;a,b 

{1, . . . ,m + l;a,b\l, . . . ,m + l;a, b)m+i-a,b 

= ^'^kiPU ■ ■ ■ ,Pm+i;Pa,Pb) +'D''''\pi, . . . ,Pm+i;Pa,Pb)- 
k^i 

where again we neglect all finite terms. The first new dipole term is given by 

1 1 



(53) 



1^'k{Pl^---^Pm+i;Pa,Pb) 



'^Pa-Pi Xik ,a 

Xm;a,6(l, ■ ■ ■ , k, . . . ,171 + l] ai,b\ ^^g^V^Jl, . . . , k, . . . , m + l] ai, b)m-a,b: 

where the Lorentz scalar a^j^ „ is given by 

Pk-Pa+Pi-Pa -Pi-Pk 



■^ik,a 



{Pk +Pi)-Pa 



(54) 



(55) 
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and the momenta of the emitter and spectator are given by 



P^i = Xik,aP'i, Pk=Pk+Pi - C^- Xik,a)p'^- (56) 

The second type of dipole on the RHS of Eq. (f53l) is given by 

V''''\pi, . . . ,Pm+i;Pa,Pb) = — 

^Pa-Pi Xi,ab 

Xm-a,b{^, . . . , mTl; ai, b\ ^^V^''% . . . , mTl; ai, b)m;a,b, (57) 

ai 

where the Lorentz scalar Xi^ab is given by 

Pa-Pb - Pi-Pa - Pi-Pb 

Xi,ab = • (58) 

Pa-Pb 

The structure of this dipole term differs from the others since it is convenient to leave the 
momentum of the spectator, pf,, unchanged. Thus the tilde momenta phase point at which 
the matrix element \1, . . . ,m + l;ai, h)m;a,b is evaluated is found as follows. The momentum 
of the emitter is parallel to pa'- 

K^ = HabP^a, (59) 

and in order to implement momentum conservation exactly all other final-state momenta 
must recoil as follows: 

where 

= P'^+Pt-Pl 

= K^+Pb■ (61) 

Note that it is all final-state particles that recoil, not just the QCD partons. Like all the 
previous dipole contributions, all particles remain on mass-shell and momentum conservation 
is exact. 



2.4.3 Summary 

We have stated that each of these dipole terms acts as a counter term for the real matrix 
element squared, |A^m+i;a,6p) in a different singular region of phase space. The sum of all 
dipole terms will therefore act as a local counter term in all singular regions of phase space. 
We need to construct a counter term for da^, where 

da"" oc d,^(™+i)(p;Q)|A^„+i;,,fc|2(p)Fj'"+^)(p), (62) 

and it is clear that such a counter term, da^, must contain a factor of the jet function that 
defines the observable, Fj. 
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Recall that the cancellation of divergences between the real and virtual contribution was 
only guaranteed for infrared safe observables defined by the following property of the jet 
function: 



F 



(m+l) 



F 



(m) 



(63) 



as one approaches an m + 1-parton configuration that is kinematically degenerate with an 
m-parton configuration. Since da^ must match the singular behaviour of da^ and in the 



singular regions of phase space we know that Fj^^^'^ 



;.(m) 



it appears we have a choice 



when constructing the counter-term da^: it can be constructed by multiplying each of the 
dipole terms by either Fj^~^^\p) or Fj^\p). Either would give a counter term that matches 
the singularities of the real cross section, but only the latter gives a finite cross section that 
can be implemented by Monte Carlo methods, as we shall now discuss. 

An obvious requirement of the counter term is that it matches the singular behaviour 
of the real term, but does not contain any additional singularities. The dipole terms are 
all proportional to the square of the Born level matrix elemenlH evaluated at the m-parton 
phase space point given by the tilde momenta {p}. Multiplying this is a factor that diverges 
in the soft and collinear regions of phase space (1/pi-Pj)- Consider some point in m + 1- 

parton phase space for which the real matrix element is finite and Fj^^^\p) is non-zero; it 
is possible that one of the phase space mappings onto m-parton phase space produces an m- 
parton configuration arbitrarily close to a singular region of phase space such that l-M^-abl"^ 
is arbitrarily large. For such a point in phase space da^ will be finite and the first choice 
that da^ is proportional to Fj^'^^\p) will be arbitrarily large, i.e. da^ will contain extra 
('spurious') singularities. 

However recall that, by definition, the Born level contribution to any observable given 
by Eq. ([3]) is finite. Thus the jet function Fj"^^ must be defined such that as the Born level 
matrix element diverges Fj^^ 0. Therefore if we make the second choice that da^ is 

proportional to Fj"^^ (p) these extra singularities will be removed. 
Finally we define the counter term as 



where 

(P-F(-))(p) 

dipoles 



da^ = #("^+i)(ij; Q) ^ iV- F(™))(p) 

dipoles 



(64) 




V,,,,{p)Ff\p) + V-^{p)Ff'\p) + V\^{p)Ff\p) 



Y,Vt{p)Ff'\p)+V^^'\p)F]j"'\p) + {a^b) 



■ (65) 



■^Actually it is the relevant colour and spin projection of the Born level matrix element squared, but that 
does not affect this argument. 
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2.5 Integrated dipole terms 

In this section we shall discuss the terms that result from the integration of the dipole terms 
over the single parton subspace leading to the singularities. We shall not give a detailed 
discussion, but will concentrate on those aspects that will be useful for the user to include a 
new process. 

In section 12.21 we outlined the basic idea that some of the divergence from dipole terms 
with initial state partons cancels the 1/e poles of the collinear subtraction term to give the 
finite convolution term in Eq. (j22p . The rest of the divergence is contained in the insertion 
operator and we stated that this will cancel the divergence of the virtual contribution such 
that the integral 

^NLO{W(^) = f [da'^ip) + da'^ip) , (66) 

J m 

is finite. After the integration over the single particle subspace no spin correlations survive, 
so the insertion operator is an operator only in colour space. The second term on the RHS is 

da^ip) Ol =m,a,b (!,••• ,m;o, 6|I({p};e)|l,... ,m; a,b)m,a,b, (67) 
where the insertion operator is given by 



(68) 



The sums runs over both initial-state and final-state partons, and the flavour-dependent 
universal singular functions V/(e) are given by 



Vg{e) = Cf 

v.(e) = ^ + (T^^-^^^^/)7 + ^^(f + (69) 

Since the singularities of this insertion operator term cancel with those of the renormalized 
one-loop matrix element, TeVJet needs only to calculate the finite parts of each. What 
we mean by the finite parts warrants a little explanation. The insertion operator can be 
re-written as 

where Q is some arbitrary scale. Using dimensional regularization the renormalized one-loop 
matrix element in the MS scheme is typically of the form 



I>'I-.-|!f7T^(^)'(4 + ^C + OW^ (71) 
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where fi is the dimensional regularization scale (which, in common with Ref. [10], we set 
equal to the renormalization scale throughout), and Q is some scale relevant to the process - 
the precise definition of Q fixes the precise form of the functions B and C. The convention 
in TeVJet is to extract a factor of 

' > (72) 



r(l -e)\Q 

from the insertion operator and the one- loop matrix element, and to calculate the finite part 
of what is left. The reference scale Q can be chosen by the user. It can be any scale as 
long as the user is consistent, and the same scale is extracted when calculating the one- loop 
matrix element. Thus for the insertion operator counter-term TeVJet needs to calculate the 
finite part of terms of the form 



where the singular functions V/(e) have the form 



Vi{e) = ^ + - + F + 0{el (74) 



and the second term of Eq. (|73|) can be expanded as a Taylor series 



log(a) log^(a) n , , 

a' = l + + + • • • • (75) 

Thus the finite part of Eq. ([75]) is given by: 

In order to put a new process into TeVJet the user must provide code that fixes Q and returns 
the finite part of Eq. (|71|) after the common factor has been removed, i.e. 

gc. (77) 



3 TeVJet Design 

TeVJet is a framework for calculating jet cross sections in NLO QCD using the dipole sub- 
traction method [10]. It has been designed to make the inclusion of a new scattering process 
as straightforward as possible. The user must supply the usual ingredients for an NLO cal- 
culation: the Born level and real emission matrix elements, and the interference between 
the Born level and one-loop amplitudes. From these TeVJet automatically calculates the 
subtraction terms and various insertion operators; and performs the integrations over phase 
space. 

TeVJet has been written in a modular way using C++; the library is made of the following 
classed, and classes that inherit from them. The main classes are: 

''Note: all of these classes are declared inside namespace TeVJet. 
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• Calculation 

• Process 

• JetFunction 

• Tools 

• PhaseSpace 

The following classes are available to the user for producing histograms: 

• Histogram 

• HistoFile 

The following classes contain the process-dependent parts of the calculation: 

• Amplitude 

• MatrixElementSquared 

• Looplntegral 

The following classes contain the process-independent parts of the calculation: 

• SubtractionTerm 

• InsertionOperator 

• SplittingKernelP 

• SplittingKernelK 

In section 14.21 we give an example of a main program. The main driving class is the 
Calculation class, which takes, on construction, a pointer to a Process and a pointer to a 
JetFunction. The Process class contains all of the process dependent parts of the calcula- 
tion, whereas the JetFunction contains the definition of the observable(s) to be calculated. 
Figure [5] shows schematically the structure of the program. In this section we shall provide 
a description of some of these classes, with particular emphasis on the methods that are 
needed by the user. In sections I3.1H3.5I we describe the main five classes listed above. In 
sections I3.6H3.7I we describe the classes used to create histograms. In sections l3.8H3.10l we 
describe the classes that store all of the process-dependent parts of the calculation. Some 
knowledge of these is required if the user wishes to include a new process in TeVJet. The user 
should not need to know the details of the four classes that contain the process-independent 
parts of the calculation, and they are not described here. 

At this point we distinguish two types of user: 

• A user who wishes to do a new calculation by writing a new Process class. They will 
need some knowledge of all of the TeVJet classes. 
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mam. cc 



Tools 



Calculation 



JetFunction* 



Process* 



std : : vector<MatrixElementSquared*>( X 2) 
std: : vector<Ainpltude*> 
std: : vector<LoopIntegral*> 
PliaseSpace*(x2) 



Subtract ionTerni*(x 4) 
Insert ionOper at or* 
SplittingKernelP* 
SplittingKernelK* 



Figure 2: Schematic diagram showing the structure of TeVJet. 



• A user who wishes to use TeVJet to calculate rates, plot distributions and so on for 
an existing process. In this case they should only need detailed knowledge of the 
JetFunction class and the two histogram classes, Histogram and HistoFile, with 
some knowledge of the other main classes. 

3.1 The Calculation class 

This is the driving class of the calculation. It collects together all of the parts of the cal- 
culation and performs the phase space integrals and convolutions with the parton distri- 
bution functions. On construction it is passed all of the process-dependent parts of the 
calculation. It also has, as private members, pointers to all of the process-independent 
parts, i.e. the four Subtract ionTerms, and the InsertionOperator, SplittingKernelP and 
SplittingKernelK. 

Calculation methods 

Calculation(TeVJet : :Process * proc, TeVJet :: JetFunction * jet); 

This is the constructor for the Calculation class. 

void InitializeO ; 
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This initializes the beam momenta, the random number generator, the Process and, if there 
are no hadrons in the initial state, the pointers to the two PhaseSpace objects. 

The following set functions should be called to set the beam energies in units of GeV and 
the random number seed before initialize otherwise default parameters will be used. The 
default beam energies depend on the number of initial state hadrons in the Process - for 
zero hadrons the beam energies are both set to for one hadron the hadron energy is 

set to 920.0 GeV and the lepton energy to 27.5 GeV; for two hadrons the energies are both 
set to 7000.0 GeV. 

void setEbeaml (double Ebeaml) ; 
void setEbeam2 (double Ebeaml); 
void setSeeddong seed) ; 

The integral of the Born level result of Eq. ()16p can be evaluated by calling the method 

void BornLoopO ; 

and the three integrals of Eq. (1170 can be evaluated by calling the methods 

void RealLoopO ; 
void VirtLoopO ; 
void ConvLoopO ; 

The following additional set functions are available to set the number of phase space points 
to choose for these four integrals. 

void setNborn(int Nborn) ; 
void setNreaKint Nreal) ; 
void setNvirt(int Nvirt) ; 
void setNconv(int Nconv) ; 

When a new process is written (see section 13. 2p the user must choose a renormalization and 
factorization scale. This must be an infrared safe function of the momenta, e.g. the total CM 
energy in the case of e^e~ collisions. The Calculation will multiply each of these scales by 
a number that can be set by the user via the following methods: 

void setRenScaleFactor (double factor); 
void setFacScaleFactor (double factor); 

Each of these numbers is set to 1.0 by default. 

void setRenSchenie(std: : string scheme); 
void setFacScheme(std: : string scheme); 

These methods are used to set the renormalization and factorization schemes. At present 
only one renormalization scheme is available, and the only possible value of scheme is 
scheme="MSbar" . Two factorization schemes are available by setting scheme="MSbar" or 
scheme="DIS" . It is important to ensure that the same factorization scheme has been used 
to extract the pdf, i.e. the user must choose a pdf that uses either the MS or the DIS 
factorization scheme, and then should set TeVJet to use the same scheme. 
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void InitPDF(std: : string pdfname, int pdfmemset) ; 

TeVJet links to LHApdf for the parton distribution functions. The user may use this method 
to choose the pdf set, and the member of that set that they wish to use. If this is not called 
the default set (cteq61) is used. 

3.2 The Process class 

The purpose of the Process class is to pass all of the process-dependent parts (e.g. matrix 
elements) to the Calculation class. It is an abstract base class, which has a number of virtual 
methods. All of the process-dependent parts of the calculation are passed to the Calculation 
class via a class that inherits from the Process class and in which these methods have 
been over-written. The process-dependent parts of the calculation are stored in the public 
members: 

std: : vector<MatrixElementSquared*> MatrixBorn; 
std: : vector<MatrixElementSquared*> MatrixReal; 
std: : vector<Amplitude*> BornAmp; 
std: : vector<LoopIntegral*> Loops; 

We define these classes later, but for now make one comment about them: A single process 
can contain several subprocesses, implemented by different members of these vectors. For 
example the process e+e" 3 jets will contain subprocesses such as e'^e" — ^ uug. 
The process class also has two pointers to PhaseSpace objects: 

PhaseSpace * PhaseBorn; 
PhaseSpace * PhaseReal; 

Process methods 

When writing a class that inherits from the Process class, the following virtual methods 
must be written. 

virtual void Initialize ()=0; 

This is called by the Calculation method Initialize () . The user should call the con- 
structors of the two PhaseSpace classes. The user should also call the constructors of all of 
the Born level Amplitude classes, the Born and real MatrixElementSquared classes and the 
Looplntegral classes. These can either be existing classes, or the user can write their own 
(see relevant sections below). 

virtual int NhadronBeaiiis()=0; 

This should be set to for processes with no hadrons in the initial state, 1 for DIS processes, 
and 2 for hadron-hadron processes. 

virtual bool DoNLO()=0; 
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This returns a flag that may be set to false if the user only wishes to calculate the Born 
level result for a process. 



virtual double getRenScale( 

std: : vector<HepLorentzVector> Momenta) =0; 
virtual double getFacScale( 

std: : vector<HepLorentzVector> Momenta) =0; 

These methods are called by the Calculation class to provide the renormalization and 
factorization scales respectively and the user should supply code to return an infrared safe 
function of the input momenta. This could be for example the total CM energy in e~^e~ 
collisions, or the virtuality of the exchanged vector boson in DIS. The first two elements of 
the input vector are the momenta of the initial state partons, followed by the final state 
partons. 

virtual double getRef Scale ( 



As explained in section 12.51 the convention in TeVJet for defining the finite parts of the 
insertion operator and the one-loop matrix element is to extract a common factor of 



and drop the divergent terms in e. The reference scale Q can be chosen by the user and this 
is what the method getRef Scale should return. It can be any scale as long as the user is 
consistent and the same scale is extracted when calculating the one-loop matrix element. 

3.3 The JetFunction class 

The purpose of the JetFunction class is to calculate the value of the jet function Fj intro- 
duced in Eq. ([3]). It is an abstract base class, which has a number of virtual methods. The 
main method is the user defined function PlotStuff , which is called several times during 
the run. TeVJet passes each phase space point with a weight to the user, and the user must 
provide the code to calculate the observable(s) of interest. For example to calculate a differ- 
ential cross section, da/dX, one should calculate X from the final state momenta and fill a 
histogram of X with the weight. 

The method PlotStuff is called several times per 'event' in RealLoop: once for each 
m + 1-body phase space point, and once for each of the different m-body phase space points 
needed for each subtraction term (see discussion in section I2.4.3|) . In general the error on an 
observable (e.g. a single bin of a histogram) can be obtained from a weighted distribution of 
events by adding their weights in quadrature. However in the subtraction algorithm there can 
be large cancellations between the real matrix element and the subtraction terms and, since 
these are correlated, one should sum them first before using them to calculate the error on an 
observable. Therefore if one is using a histogramming package that automatically calculates 
the errors on each bin using the squares of the weights one will overestimate the errors if the 



std: : vector<HepLorentzVector> Momenta) =0; 




(78) 



24 



real matrix element and each subtraction term are filled separately. In order to avoid this 
there is a method Terminate, which is called at the end of each full event. The user should 
fill a temporary histogram in PlotStuf f each bin of which is added to the main histogram in 
Terminate. The same applies to the convolution term calculated by ConvLoop - PlotStuf f 
is called twice with two correlated phase space points followed by one call to Terminate. 

Since the details of the observable (s) to be calculated arc contained in a class that inherits 
from the JetFunction class and in which these methods have been over- written, we describe 
the methods below. 

JetFunction methods 

The constructor of the JetFunction class takes a std: : string, which is the name of the 
output file: 

JetFunction(std: : string filename) ; 
The JetFunction class has a standard vector of pointers to Histogram objects: 

std: : vector<Histogram*> m_histovec; 
and a pointer to a HistoFile object: 

HistoFile * m_outfile; 

The following method will write all of the histograms stored in mJiistovec to the output 
file: 

void WriteFileO; 

When writing a class that inherits from the JetFimction class, the following virtual methods 
must be written. 

virtual bool KeepPhasePoint ( 

std: : vector<HepLorentzVector> finals)=0; 

This method is called after the generation of each phase space point, before the weight is 

calculated. If this returns true the point is kept and the weight calculated, if it returns 
false the point is rejected and the weight never calculated. If the user wishes to put some 
cuts on phase space they should be implemented here. 

virtual void PlotStuf f (double weight, 

std: : vector<HepLorentzVector> finals, 
std: : string type)=0; 

This method is called several times by the Calculation class. The variable type indicates 
which type of contribution is being calculated. It is called once for each phase space point 
in BornLoop and VirtLoop (type is born or virtual respectively); several times for each 
phase space point in RealLoop - once for the real emission ('type is real) and once 
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for each dipole configuration (type is subtraction); and twice for each phase space point 
in ConvLoop (type is convolution). Here the observable(s) should be calculated using the 
final state momenta and the weight. For example to calcTilatc a differential cross section, 
da/dX, one should calculate X from the final state momenta and fill a histogram of X with 
the weight. 

virtual void Terminate ()=0; 

As explained above, the phase space points for the real emission \M.\'^ and the dipole con- 
figurations are correlated and one will over-estimate the error if one fills the histogram once 
for this 'event' and once for each 'counter-event'. One way to avoid this is to fill a dummy 
histogram in PlotStuf f in this way, and at the end of the 'full event' to copy this histogram 
into the main histogram. The method Terminate () gets called once per real phase space 
point in RealLoop for this purpose. It is also called once per full 'event' in ConvLoop. 

For this purpose the following method exists to copy the contents of one histogram 
(histo2) into another histogram (histol): 

void FillEvent (Histogram * histol. Histogram * histo2) ; 
3.4 The Tools class 

All parts of the program need access to the Standard Model parameters (for which, we note, 
we use natural units based on GeV, with c = 1). We therefore provide them through the 
Tools class, which is a 'singleton' class, as we now describe. 

The singleton pattern 

The 'singleton pattern' is essentially the object oriented way of declaring a global variable. 
A singleton class has a private constructor, but also has as a static private member a pointer 
to itself that can be accessed by a public method. Here is a skeleton of the class declaration: 

class Tools{ 
public : 
"Tools ; 

static Tools * getlnstanceO ; 

private : 

//private constructor 
Tools ; 

static bool instance; 
static Tools * g_tools; 

};// end class Tools 

The method getlnstance returns the private pointer (on the first call the constructor is 
called): 
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TeVJet: :Tools * TeVJet: : Tools : :getlnstance(){ 
if ( ! instance) { 

g_tools = new TeVJet: : Tools () ; 
instance = true; 
return g_tools; 

} 

else{ 

return g_tools; 

} 



The Tools class is therefore a safe global variable that stores the values of the SM parameters. 
It does not get created until one needs it for the first time, at which point the constructor is 
called and the default values for the parameters are set. 



Tools methods 

Through the methods of the Tools class the user may access the parameters, or set them to 
values other than the defaults. 

double Alpha_s(int Nloop, double scale); 

This returns the value of as{l-t), evaluated at renormalization scale scale, where the running 
uses the /3-function evaluated at either one- or two-loop order (Nloop=l or 2 respectively). 

double Alpha_ em (double scale) ; 

This returns the value of oemC/^), evaluated at renormalization scale scale. The running 
of OEM with scale is not implemented at present, this method simply returns the value of 
the coupling at a renormalization scale equal to the Z boson mass, i.e. aEM(-^z)- However 
since this coupling is always accessed through this method, the running of the coupling could 
easily be included if needed. 

The following 'get' methods may be used to obtain the values of various parameters. 

double getUpMassO; 

This returns the up quark mass. 

double getDownMass ; 
This returns the down quark mass. 

double getCharmMassO ; 
This returns the charm quark mass. 

double getStrangeMass ; 
This returns the strange quark mass. 
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double getTopMass ; 
This returns the top quark mass. 

double getBottomMass () ; 
This returns the bottom quark mass. 

double getTopWidthO ; 
This returns the top quark width. 

double getBottomWidthO ; 
This returns the bottom quark width. 

int getNoOf LightFlavours ; 
This returns the number of hght quark flavours. 

double getWmassO; 
This returns the W boson mass. 

double getZmass ; 
This returns the Z boson mass. 

double getHmassO; 
This returns the Higgs boson mass. 

double getWwidthO ; 
This returns the W boson width. 

double getZwidthO; 
This returns the Z boson width. 

double getHwidthO ; 
This returns the Higgs boson width. 

double getSin2W() ; 
This returns the value of sin^ 9w- 

double getGfermiO; 
This returns the value of Fermi's constant, Gp- 

double getAlphaMZO ; 

28 



This returns the value of the electromagnetic coupling at the Z boson mass, a^yiiMz)- 

double getAlphasMZO ; 
This returns the value of the strong coupling at the Z boson mass, as{Mz)- 

int getNcO; 
This returns the number of colour states of quarks, N^. 

double getCFO ; 
This returns the value oiCp. 

double getCAO ; 
This returns the value of Ca- 

double getTRO ; 

This returns the value of Tr. 

The parameters used by TeVJet can be altered by the user, note however that these 
parameters cannot be changed during the run, and should be set in the main program before 
the constructors for any of the TeVJet classes are called. This is because these constructors 
may need access to the SM parameters. For this reason there is a method 

void FreezeO ; 

which 'freezes' the parameters stored in the Tools class so that subsequent calls to any of 
the 'set' methods will not alter the parameters and will print a warning message. This is 
called by the JetFunction and Process constructors. 

The default value for each of the parameters is listed in Table [H The following 'set' 
methods may be used to set the values of these parameters. 

void setUpMass (double mass); 

Sets the up quark mass to mass. 

void setDownMass (double mass); 
Sets the down quark mass to mass. 

void setCharmMass (double mass); 
Sets the charm quark mass to mass. 

void setStrangeMass (double mass); 
Sets the strange quark mass to mass. 

void setTopMass (double mass); 
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Parameter 


Name 


Type 


Default 


ruu 


m_UpMass 


double 


0.0 


GcV 


rrid 


ni_DownMass 


double 


0.0 


GeV 


rUc 


m_CharniMass 


double 


1.25 


GeV 


rus 


m_StraiigeMass 


double 


0.0 


GeV 


mt 


m_TopMass 


double 


174.3 


GeV 


mi, 


ni_BottomMass 


double 


4.7 


GeV 


Tt 


m_TopWidth 


double 


1.508 


GeV 


Tb 


m_BottomWidth 


double 


0.0 


GeV 


Nf 


m_NoOf LightFlavours 


int 


5 




EW scheme 


ni_E W c n V e nt i on 


int 


3 




/ 71 /T \ 


m_alpha 


double 


1/128.89 




Gf 


m_Gf ermi 


double 


1.16639 X 10~^ 


GeV ^ 


Sin d-w 


m_sin2thetaW 


double 


0.231 




Mw 


m.Wmass 


double 


80.41 


GeV 


Mz 


m_Ziiiass 


double 


91.188 


GeV 


Mh 


m_Hmass 


double 


120.0 


GeV 


i w 


m_wwiatn 


double 


2.124 


GeV 


Tz 


m.Zwidth 


double 


2.4952 


GcV 


m-r 


m_TauMass 


double 


1.777 


GeV 


as{Mz) 


m_alpha_s_MZ 


double 


0.118 




Nc 


m_Nc 


int 


3 




Tr 


m_TR 


double 


0.5 





Table 1: Default values for the SM parameters in TeVJet. 



Sets the top quark mass to mass. 

void setBottomMass (double mass); 
Sets the bottom quark mass to mass. 

void setTopWidth (double width) ; 
Sets the top quark width to width. 

void setBottomWidth (double width); 
Sets the bottom quark width to width. 

void setNoOf LightFlavours (int N) ; 

Sets the number of light quark flavours to N. 

The electroweak parameters require a little explanation. There are 6 parameters, which 

can be chosen to be the SU(2) coupling strength g, the weak mixing angle sin^ 6y/, the 
electromagnetic coupling constant qem and the masses of the W , Z and H bosons. However 
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m_EWconvention 


Input parameters 





OEuiMz), Gf, sin^6'iy 


1 


Mw, Gp, sin^ 6w 


2 


Mz, aEM(Mz), sin^Ow 


3 


Mz, Mw, Gf 



Table 2: The input parameters that can be set by the user for different m_EWconvention 
values. The defaults for these parameters are listed in Table [H Note the default values are 
only used for the input parameters and the remaining parameters are calculated from these. 



tree-level gauge invariance imposes two constraints such that there are only 4 free parameters, 
and as such one should not allow all six parameters to be set independently. TeVJet follows 
the convention in ALPGEN [2], which does the following. It treats the Higgs mass as a free 
parameter, it then takes 3 input values, which can be set to the best fit to experimental data 
and calculates all other parameters from these using the constraints from tree-level gauge 
invariance. There is a choice of which 3 parameters are used as the input, which is controlled 
by an integer variable m_EWconvention and can be set via the method 

void setEWconvention(int i) ; 

This is equivalent to the variable iewopt in [2]; the possible values are 0, 1, 2 or 3. Table [2] 
lists the input parameters for the different values of ni_EWcorLveiitioii. 
The Higgs boson mass can be set via the method 

void setHmass (double mass); 

which sets the Higgs boson mass to mass. 

The following methods can be used to set the other 3 input parameters. 

void setWmass (double mass); 
Sets the W boson mass to mass. 

void setZmass (double mass); 
Sets the Z boson mass to mass. 

void setSin2W(double sin2w) ; 
Sets the value of sin^ 6w to sin2w. 

void setGfermi (double Gfermi) ; 
Sets the value of Fermi's constant to Gfermi. 

void setAlphaMZ (double alpha); 
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Sets the value of the electromagnetic coupling at the Z boson mass, a^y[{Mz)- 

Once one of the three input vahics for the current m_EWconvention value is set, the other 
parameters are automatically recalculated. If one attempts to set one of the parameters that 
is not an input parameter for the current in_EWconvention value, TeVJet will print a warning 
statement, and that parameter will not be changed. 

void setWwidth (double width) ; 
Sets the W boson width to width. 

void setZwidth (double width) ; 
Sets the Z boson width to width. 

void setHwidth (double width) ; 
Sets the Higgs boson width to width. 

void set AlphasMZ (double alphas); 
Sets the value of the strong coupling at the Z boson mass, asiMz)- 

void setNc(int N) ; 

Sets the number of colour states of a quark to N. If this is set the values returned by the meth- 
ods getCFO and getCAO are altered accordingly. Also the QCD /3-function is recalculated 
and the running of as will change. 

void setTR (double TR) ; 

Sets the value of Tr to TR. 

3.5 The PhaseSpace class 

There are currently three phase space classes that inherit from the PhaseSpace class: 
TwoBodyPhaseSpace, ThreeBodyPhaseSpace and FourBodyPhaseSpace. 

An element of massless two-body phase space can be written as 



where Q* and are the polar and azimuthal angles of one of the particles in the centre of 
mass frame. By default, points in two-body phase space are chosen flat in cos0* and For 
processes with a hadron in the initial state the partonic cross section must be folded with 
the parton distribution functions, so we need to calculate quantities like 



d(j)^'^\pi,P2;Q) = — d cos 9* d(l), 



(79) 



I 




(80) 
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In general it is very inefficient to pick flat in r]. Instead, we can use the mapping 

drj d cos 6* = — —dptdy, (81) 
Vs 

where pt and y are the transverse momentum and rapidity of one of the particles in the centre 
of mass frame of the beams. It is possible to pick flat in y and as via the method SetMode. 

The three-body phase space generator uses a two-body configuration and the phase space 
factorization of Eq. ([9]) to produce one additional particle. In 4 dimensions the dipole phase 
space factorization is 

d(t>^^Hpi,Pj,Pk;Q) = ^^Q^3 (1 - 0(2j(l - 6(7/ij,fc(l -yij,k))dzidyij^kd(l) 

xd4>^^\p^j,Pk;Q). (82) 

Since this is an exact factorization, in order to integrate a well-behaved function of phase 
space it would be sufficient to follow the following algorithm to produce three-particle final 
states: 

1. Pick a two-body configuration {pi2,'Pz) with weight W2 

2. Pick the dipole variables ij, yij^k and (j) in the appropriate ranges 

3. Use these and the two-body configuration to reconstruct the three-body configuration 
{PiiP2-,Pci) where pi2 is the emitter and p^ the spectator. This should be given weight 

-^3 = -W2 X ^^^3 (1 - j/12,3)- (83) 

However the integration of the real matrix element minus the subtraction terms over phase 
space typically contains square root singularities in the dipole variables, i.e. 

\M^^i?Ff'^'\) - Y: {T^-F!r%)] m,,k - (84) 



dipoles 



and 



dipoles / L 

Therefore in order for the Monte Carlo integration to converge one must introduce a mapping 
to remove these singularities. In order to do this we introduce a Jacobian factor of the form: 



J= ^ r- (86) 

i.e. the dipole variables are chosen such that the weight should be multiplied by this factor, 
which tends to zero as any of the limits yij k — > or 2j — > 0, 1 are approached. Since the 
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integrand contains square root singularities in these limits for any choice of i,j and k TeVJet 
uses a multichannel phase space generation so that instead of fixing the choice of emitter and 
spectator all possible choices are used. Each channel is chosen with equal probability, and 
the total Jacobian factor J' is given by 



where J7i is the Jacobian factor for the ith. channel, given by Eq. (I86p . 

Four-body configurations are generated from three-body configurations in a similar way. 

PhaseSpace methods 

The constructor for the PhaseSpace classes takes a std : : vector with the masses of the final 
state particles. 

PhaseSpace (std: : vector<double> massvec) ; 

At present only massless phase space generation is implemented, so massvec must be filled 
with zeros, but we anticipate the inclusion of masses soon. The way that two-body phase 
space is picked can be altered via the method: 

void SetMode(std: : string mode); 

Possible values of mode are lepton-lepton where two-body phase space is picked fiat in 
9* and (j); DIS where two-body phase space is picked fiat in the rapidity y and as some 
power of the pt of one of the outgoing particles in the centre of mass frame of the 

beams; or hadron-hadron where two-body phase space is picked fiat in the rapidities of the 
two outgoing particles, yi^2, and the transverse momentum of the particles to some power, 
^ptpow^ in the centre of mass frame of the beams. 

When the mode has been set to DIS the user may choose the minimum pt and the power 
ptpow via the methods 

void SetPtpow (double ptpow) ; 
void SetPtmin (double ptmin) ; 

3.6 The Histogram class 

The Histogram class can be used to store the distributions calculated by the JetFunction 
class. Two versions of the Histogram class source code are available: HistogramSimple.ee 
and HistogramRoot.ee. The 'simple' version is a very basic histogram package and is used 
by default, the other version is a wrapper for the Root [16] class THIF. In order to use the 
Root version one must remove HistogramSimple.ee from the source directory and replace 
it by the Root version, which is in the directory Hist code/; uncomment the relevant lines 
in the header file (Histogram. h); and uncomment the lines in the Makefile that link to the 
Root libraries. If the user wishes to use another histogram package they should write the 
necessary source code for the methods described below. 




(87) 



34 



Histogram methods 

The constructor takes the name and title of the histogram, the number of bins, and the 
minimum and maximum values of the quantity: 

HistogramCstd: : string name, std::string title, 
int Nbins, double min, double max); 

The following methods are available to fill quantity x either unweighted (i.e. with weight 1) 
or with weight weight: 

void Fill (double x) ; 

void Fill (double x, double weight); 

If the following method is called, the Histogram will store the error on each bin by adding 

the weights in quadrature: 

void Sumw2() ; 

The following methods return the width, error, centre or content of bin number bin (Note: 
the bins are labeled from 1 to Nbins): 

double GetBinWidth(int bin) ; 
double GetBinError (int bin) ; 
double GetBinCenter (int bin) ; 
double GetBinContent (int bin); 

The following method returns the number of bins: 

int GetNbinsXO; 

The following methods return the minimum or maximum values of the observable: 

double GetMinimumO ; 
double GetMaximumO ; 

The following method returns the title of the histogram: 

std::string GetTitleO; 

3.7 The HistoFile class 

The HistoFile class is the file that the Histograms stored in the vector mJiistovec of 
the JetFunction class are written to at the end of the run. Two versions of the HistoFile 
class source code are available: HistoFileSimple.ee and HistoFileRoot.ee. The 'simple' 
version is a text file that contains names of each histogram and the values and errors of 
each bin and is used by default. The other version is a wrapper for the Root [16] class 
TFile. In order to use the Root version one must remove HistoFileSimple . cc from the 
source directory and replace it by the Root version, which is in the directory Histeode/; 
uncomment the relevant lines in the header file (HistoFile .h); and uncomment the lines in 
the Makefile that link to the Root libraries. If the user wishes to use another file format they 
should write the necessary source code for the methods described below. 
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HistoFile methods 



The constructor takes a std: : string, which is the name of the output file. This will be 
appended with .simple or .root depending on which version of HistoFile. cc is used: 

HistoFile (std :: string fname) ; 
The following method writes the Histograms in histos to the output file: 

void Write (std : :vector<Histograiii*> histos); 
The following method closes the output file: 

void CloseO ; 

3.8 The MatrixElementSqueired class 

The MatrixElementSquared class is an abstract base class. The mechanism for providing 
code to calculate the square of the matrix element for a sub-process that contributes to a 
process is to make a class that inherits from the MatrixElementSquared class. Pointers to 
instances of this then get 'pushed back' into one of the vector members of the Process in its 
Initialize method. When writing a class that inherits from the MatrixElementSquared 
class, the user must initialize the following members: 

int m_NInitColour ; 

This is the number of initial state colour states, e.g. for the process qg qg this would be 
3 X 8 = 24. 

int m_BoseSymFactor ; 

This is the Bosc symmetry factor that one must divide by when there are identical particles 
in the final state. 

In addition to these, the following members must be initialized: 

std: :vector<int> m_InitialID; 
std: :vector<int> m_FinalID; 
std: :vector<int> m_IDlist; 

These are initialized in the methods SetParticlelDs and FilllDlist (see below). 
MatrixElementSquared methods 

When writing a class that inherits from the MatrixElementSquared class, the following 
virtual methods must be written. 

virtual void SetParticlelDs ()=0; 
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In this method the user must push back the PDG ID codes of the two initial state particles 
into the member object ni_InitialID and the ID codes of all final state particles into the 
member object m_FinalID. Finally the user should call the method FilllDlist. Note: this 
method must be called by the constructor of this class. 

virtual double GetMSquared( 

HepLorentzVector parti, 
HepLorentzVector part2, 

std: : vector<HepLorentzVector> finalstate, 
double RenScale)=0; 

This is a method that should return the square of the matrix element given the initial and 
final state momenta, and the rcnormalization scale. The momenta will be passed in the order 
set by the user in the method SetParticlelDs. 

3.9 The Amplitude cleiss 

The Amplitude class is an abstract base class, which has a number of members, which must 
be initialized and virtual methods, which must be overwritten. 

It has been shown [17] that QCD amplitudes can be decomposed (the colour-flow decom- 
position) in terms of different colour structures as 

M=J2^aAa, (88) 
a 

where the factor associated with each colour structure, Aa, is called a partial amplitude. 
Thus the square of the matrix element for a scattering process can be written as 

\Mf = J^CaAaClAl 

a,b 

= ^C-^bgafe, (89) 

a,b 

where the elements of the colour flow matrix are given by Cab = Cacl and the 'quark matrix' 
is given by Qab = AaAl. Many of the terms needed in the dipole subtraction method are 
proportional to the square of the colour correlated Born level matrix element given by a 
similar expression, but with different colour matrices: 

\M''\' = T.^abQab- (90) 
a,b 

Thus the Calculation class needs access to both 

• the 'quark matrix', Q^b, and 

• the set of colour matrices, {C**-^}. 
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In fact since we are usually interested in summing/averaging over final-/initial-state spins, 
the Calculation class requires the object 

Qab = T.Qab, (91) 

where the sum is over the set of helicity configurations of the external particles. In addition 
to this, in the case of a gluon emitter, the subtraction terms have spin correlations. For 
scattering processes with an external gluon, each partial amplitude can be written as: 

•^a, — ^ix'^ai (92) 

where e* is the polarization vector of the external gluon. The spin averaged quark matrix is 
then given by: 

Qah = I^A^b 

{s} 

= E E ^T^i'A>:.Ai\ (93) 

{S'} Sg = l,2 

where in the last line the first sum is over the set of helicity configurations of all external 
particles except the gluon. To obtain the spin averaged quark matrix one would do these 
sums, the second of which can be done by the replacement 

E ^T^l' ^ -5m- (94) 

Sg = l,2 

In the case of a gluon emitter, the subtraction term is found by the replacement 

^T^y ^ V- (95) 

Thus the Calculation class requires access to: 

• Each partial amplitude projected onto the helicity space of each external gluon, Aa- 

When writing a class that inherits from the amplitude class, the user must initialize the 
following members: 

std: :vector<int> m_InitialID; 
std: :vector<int> ni_FinalID; 
std: :vector<int> m_IDlist; 

These are initialized in the methods SetParticlelDs and FilllDlist (see below). 

unsigned int m_Nhel; 
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When calculating the subtraction term with a gluon emitter, the Calculation class explicitly 
performs the sum/average over final/initial state spins. To calculate the subtraction term 
with a particular gluon emitter, it needs to know the total number of helicity configurations 
of the other external particles, this is m_Nhel. 

std: :vector<std: :vector<int> > m_HelicityVec ; 

It also needs to know each of these helicity configurations. This object stores all of these. 

unsigned int m_NinStates ; 

This is the total number of initial state colour and spin states for the sub-process. For 
example for the process qg qg there are 2x2 = 4 spin states and 3 x 8 = 24 colour states, 
so m_NinStates should be set to 4 x 24 = 96. 

unsigned int m_Ncolour; 

This is the number of terms in the decomposition of Eq. ([88]), i.e. the number of different 
colour structures. 

HepMatrix m_matrix; 

This is the 'colour flow matrix' of Eq. ()89p . 

Amplitude methods 

When writing a class that inherits from the Amplitude class, the following virtual methods 
must be over-written. 

The constructor of the class that inherits from the Amplitude class must call the first 
three methods: SetParticlelDs, SetHelicityConf igs and InitColourMatrices. 

virtual void SetParticlelDs () ; 

In this method the user must push back the PDG ID codes of the two initial state particles 
into the member object m_InitialID and the ID codes of all final state particles into the 
member object m_FinalID. Finally the user should call the method FilllDlist. 

virtual void SetHelicityConf igs () ; 

When calculating the subtraction term with a gluon emitter, the Calculation class explicitly 
performs the sum/average over final/initial state spins. To calculate the subtraction term 
with a particular gluon emitter, it needs access to each partial amplitude projected onto the 
helicity space of that gluon for each helicity configuration of the other external particles. 
The two helicity states of quarks and gluons are both stored as ±1. In this method, the user 
should 'push back' into the vector m_HelicityVec a set of vectors of integers (either ±1), 
which represent the helicity configurations of the other external particles. 

virtual void InitColourMatrices () ; 
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Here the user should mitialize the elements of the colour flow matrix, m_matrix. The user 
should also initialize the set of colour flow matrixes, {C'*'^}, which are needed for the sub- 
traction terms. 



virtual HepMatrix GetAmplitudeMatrixQuark( 

std: : vector<HepLorentzVector> momenta, 
double RenScale)=0; 

This returns the 'quark matrix', Qab, of Eq. ([89]) . 

virtual complex<HepLorentzVector> GetAmplitudeGluon( 



This should return the Ncol-th partial amplitude projected onto the helicity space of the 
gluon, which is the ijtida-th element of m_IDlist. The other external helicities are given 
by the ihel-th element of m_HelicityVec. 

virtual HepMatrix GetColourMatrix(int ijtilda, int ktilda)=0; 



3.10 The Looplntegral class 

The Looplntegral class is an abstract base class. The mechanism for providing code to 
calculate the interference between the Born level and 1-loop amplitudes for a sub-process that 
contributes to a process is to make a class that inherits from the Looplntegral class. Pointers 
to instances of all such classes then get 'pushed back' into the vector in the Initialize 
method of the Process class. 

Looplntegral methods 

The Looplntegral class has one virtual method, which must be over-written: 
virtual double GetLoop(HepLorentzVector parti, 



The user should provide code that returns the finite part of the renormalized one-loop matrix 
element. This should be calculated in the MS scheme. What we mean by 'the finite part' 
of the renormalized one-loop matrix element is explained in section 12.51 Using dimensional 
regularization the one-loop matrix element in the MS can be written in the form 



std: : vector<HepLorentzVector> momenta, 
double RenScale, 

int ijtilda, int Ncol, int ihel)=0; 



This returns the colour flow matrix C*-''^. 



HepLorentzVector part2, 

std: : vector<HepLorentzVector> finalstate, 
double RenScale) =0; 



l-^l?-ioop 



as 1 f 47r/u 



27r r(l - e) V 




{ 



(96) 
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where fi is the renormahzation scale, and Q is some scale relevant to the process. The user 
should provide code that returns 

—C. (97) 

The scale Q must be the scale that is returned by the getRef Scale method of the Process 
class (see section 13. 2p . 

4 Using TeVJet 

4.1 Installing TeVJet 

The TeVJet source code is available from: 

http : / / www . hep . man . ac . uk/u/chris/tevj et 
To unpack this type: 

tar -zxvf TeVJet-1 .O.O.tar .gz 

This will create a directory TeVJet-1 . . 0. TeVJet must be linked at compilation time to 
HELAS [18], CLHEP [19] (version 1.8 only at present) and LHAPDF [20]. The source code 
for the HELAS library is distributed with TeVJet, and can be compiled by typing (in the 
directory TeVJet -1 .0.0) 

make libdhelasS 

CLHEP (version 1.8) and LHAPDF must both be installecll and the user must edit TeVJet's 
Makefile to provide the path to the libraries. Once this is done type: 

make 

This will build the TeVJet library and some example programs. For processes with initial 
state hadrons TeVJet needs access to the LHAPDF grid files type: 

make PDFsets 

This will create a symbolic link to the grid files. 

One of the example programs that will have been built is called EEtoSjet. We show here 
its main program and the output that it produces. 

4.2 Example main program 

//standard lib headers 
#include <iostream> 
#include <cmath> 

//my lib headers that are always needed 

^They can be obtained from the web addresses given in Refs. [19,20] together with installation instructions. 
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#include "Tools. h" 
#include "Calculation. h" 
#include " JetFunction.h" 
#iiiclude "Process. h" 

//process specific user headers 

#include "ProcEEtoSjet .h" 
#iiiclude "UserEEtoSjet .h" 

int main(int argc, char** argv) { 

long seed = 12345; 

std: : string f ilename="EEto3jet_12345" ; 



//set the SM parameters 

TeVJet: : Tools: :getInstance()->setEWconvention(3) ; 

TeVJet: : Tools: :getInstance()->setWmass(80.41) ; 

TeVJet: : Tools: :getInstance()->setZmass (91 . 188) ; 

TeVJet : : Tools : : getlnstance () ->setGf ermi (1 . 16639*pow( 10 , -5) ) ; 

TeVJet: :Tools: :getInstance()->setAlphasMZ(0. 118) ; 



//set some parameters 
int Nborn = 1000 

int Nvirt = 1000 

int Nreal = 1000 



//no of 'events' 



double CMenergy = TeVJet : :Tools : :getlnstance()->getZmass() ; //GeV 
bool useZ = false; 

double Ebeaml = CMenergy/2; 
double Ebeam2 = CMenergy/2; 

double f ren = 1.0; 

TeVJet :: Process * proc; 
TeVJet :: JetFunction * jet; 

//process specific stuff 

proc = new TeVJet :: ProcEEto3 jet (CMenergy ,useZ) ; 

jet = new TeVJet : :UserEEto3 jet (filename .CMenergy ,useZ) ; 

//process specific ends 

TeVJet: : Calculation calc(proc, jet) ; 
//user can set stuff here 
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calc . setNborn(Nborn) ; 
calc . setNreal (Nreal) ; 
calc . setNvirt (Nvirt) ; 
calc . setSeed(seed) ; 
calc . setEbeaml (Ebeaml) ; 
calc . setEbeam2 (Ebeam2) ; 

calc . setRenScaleFactor (f ren) ; 

//initialize stuff 
calc. Initialize ; 

//calculate born level 
calc . BornLoop ( ) ; 

//if required calculate NLO contribution 
if (calc . m_proc->DoNLO ) { 

calc .RealLoopO ; 

calc.VirtLoopO ; 

if (calc .m_proc->HadronBeams ){ 
calc . ConvLoop ( ) ; 

} 

} 

delete jet; 
jet=0; 

delete proc; 
proc=0; 

return 0; 



4.3 Example output 



TeVJet Version 1.0.0 - Written by C.Tevlin 
A general framework for calculating jet cross 
sections in NLO QCD. This is an implementaion 
of the Dipole Subtraction Method of S.Catani & 
M . H . Seymour . 
(Last Updated 1/2/2008) 
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Info: TeVJet: : Tools: :PrintSMPaxameters() 

Table of SM parameters used for this run: 



UpMass 

DownMass 

CharmMass 

StrangeMass 

TopMass 

BottomMass 

TopWidth 

BottomWidth 

No of light flavours 

EW scheme 

alpha_ew(MZ) 

Fermi's constant 

sin"2(thetaW) 

Wmass 

Zmass 

Hmass 

Wwidth 

Zwidth 

TauMass 

alpha_s(MZ) 

Nc 

TR 






1.25 



174.3 

4.7 

1.508 



5 

3 

0.00755099 

1.16639e-05 

0.222421 

80.41 

91.188 

120 

2.124 

2.4952 

1.777 

0.118 

3 

0.5 



Info: TeVJet: :Calculation: : InitializeO 

Additional parameters for this run: 



Ebeaml = 45.594 

Ebeam2 = 45.594 

Seed = 12345 

Renorm Scheme = MSbar 

Renorm Scale Factor = 1 



Info: Calculating the Born level 
i = 
i = 1000 
i = 2000 
i = 3000 
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i = 4000 

i = 5000 

i = 6000 

i = 7000 

i = 8000 

i = 9000 

Info: Calculating the Real Diagrams and Subtraction Terms ... 

i = 

i = 1000 

i = 2000 

i = 3000 

i = 4000 

i = 5000 

i = 6000 

i = 7000 

i = 8000 

i = 9000 

Info: Calculating the Virtual + Insertion Operator ... 

i = 

i = 1000 

i = 2000 

i = 3000 

i = 4000 

i = 5000 

i = 6000 

i = 7000 

i = 8000 

i = 9000 



Average value of (1-T) NLO 

order alpha_s contribution 
order alpha_s'"2 contribution 



= 0.0558852 +- 0.00224712 

= 0.338854 +- 0.00565262 
= 1.14195 +- 0.154111 



By default the only quantity calculated is the average value of 1-thrust, which can be calcu- 
lated as a perturbative series in a^: 

(1 - T) = aoa, + aia'i + ... (98) 
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This program prints the full NLO QCD prediction and the two coefficients ao and ai. We 
note that this example program is a fairly low statistics run intended for validation purposes, 
and we have also calculated these quantities with much higher statistics, finding agreement 
with previous calculations (see for example [5]). 

The program also produces an output file EEto3jet_12345 . simple, which contains the 
histograms. This is a text file that is too long to show here, but we note that our version 
may be found in the directory ExampleOutput/. 



5 Results 



5.1 



e+e' 



3 jets 



In this section we present some results for three-jet event shapes at LEP. We have used the 
one-loop renormalized matrix element from [7] . The renormalization scale is set to the centre 
of mass energy, which is chosen to be equal to the Z boson mass (which is set to 91.188 GeV). 
We present results for the thrust distribution and the C-parameter [4]. 

The differential cross section for thrust can be written, at next to leading order, as 



-(1-t 



da 



27r 



27r 



At{t)2TTbQ log 



+ Bt{t) 



(99) 



where ^ is the renormalization scale, s is the square of the centre of mass energy, ctq is the 
leading order cross section for e^e~ — > hadrons, and 6o is given by 



33 - 2N 



f 



12tt 



(100) 



The two coefficient functions At and Bt are universal and perturbatively calculable. Figure 
[3] shows the two coefficient functions, and Figure [H shows the full next to leading order 
QCD prediction (we have used as{Mz) = 0.118). The differential cross section for the C- 
parameter can be written in a similar form in terms of two coefficient functions Ac and Be- 
Figure [5] shows these two coefficient functions, and Figure [6] shows the full next to leading 
order QCD prediction for the C-parameter. 



5.2 ep — ^ e + 1 jet 

We have also implemented the process ep — > e + 1 jet. The beam energies chosen were a 
50 GeV electron beam on a 500 GeV proton beam (i.e. s = 10^ GeV). We have used the 
parton distribution function CTEQ5M [21] and the renormalization and factorization scales 
are set to Q where = —q^ and q is the momentum transfer between the electron and 
the proton. A cut on the transverse momentum of the electron {pt > 10 GeV) was applied 
and to define the jets we have used the kt algorithm [22] as implemented in KtJet [23] 
in its 'inclusive' mode [24] with R = 1. Figure [7] shows the transverse momentum and 
pseudorapidity distributions of the electron and the leading jet in the electron-proton centre 
of mass frame. Good agreement with the results of Ref. [25] is noted. 
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Figure 3: The distributions for the thrust coefficient functions (a) At and (b) Bf. 
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Figure 4: The full next to leading order QCD prediction for the thrust distribution. 
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Figure 5: The distributions for the thrust coefficient functions (a) Ac and (b) Be- 
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Figure 6: The full next to leading order QCD prediction for the C-parameter. 
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Figure 7: The transverse momentum and pseudorapidity distributions of the electron and 
leading jet in the electron-proton centre of mass frame. A cut on the transverse momentum 
of the electron {pt > 10 GeV) has been applied, and for the pseudorapidity distribution of 
the leading jet a cut on the pt of the jet pt > 15 GeV was applied. 
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6 Conclusion and Future Plans 



We have presented the Monte Carlo program TeVJet for calculating jet observables in NLO 
QCD, which is a direct implementation of the dipole subtraction method [10]. In order 
to implement a new scattering process the user must provide code to calculate the matrix 
element squared for the Born level and real emission diagrams, and the interference between 
the Born and 1-loop amplitudes. In fact the user must provide the colour algebra necessary 
to evaluate the square of the colour correlated Born level matrix element, as well as the Born 
level amplitude projected onto the helicity space of each external gluon. TeVJet automatically 
calculates the subtraction terms and the operators I, P and K, which are required by the 
dipole subtraction method; and performs the multi-channel phase space integrations and the 
convolution with the parton distribution functions in the case of initial state hadrons. 

We have demonstrated with the example processes e~^e~ — > 3 jets and ep ^ e + 1 jet that 
the framework works and that the results agree with existing calculations. 

Our immediate plans are to extend the Monte Carlo program to include the possibility 
of massive external partons. Finally we note that it should be possible to fully automate 
the generation of code to calculate the Born level and real matrix elements so that in future 
the only process dependent information the user must provide is the virtual term, i.e. the 
interference between the Born and 1-loop amplitudes. 
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